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Abstract. We propose a numerical method to approximate the solution of 
second order elliptic problems in nonvariational form. The method is of 
Galcrkin type using conforming finite elements and applied directly to the 
nonvariational (nondivergence) form of a second order linear elliptic problem. 
The key tools are an appropriate concept of "finite element Hessian" and a 
Schur complement approach to solving the resulting linear algebra problem. 
The method is illustrated with computational experiments on three linear and 
one quasilinear PDE, all in nonvariational form. 



1. Introduction 

Finite element methods (FEM) arguably constitute one of the most successful 
method families in numerically approximating elliptic partial differential equations 
(PDE's) that are given in variational (also known as divergence) form. 

For the reader's appreciation of this statement we briefly introduce standard 
FEM concepts. Let Q be a given domain (open and bounded set) in R d , d G N, 
/i o-a,/3 — o-i3,a : — R, be given functions with the appropriate regularity such that 
the operator div(AVu), for A :— [a a ,p\a,p=x,...,d-i makes sense, is elliptic and there 
is a unique function u : fi —> R satisfying div (AVu) = f with u = 0on dQ [GT83, 
for details]. The classical solution, u, of this problem can be characterized by first 
writing the PDE in weak (also known as variational) form using Green's formula: 

u G & and satisfies a{u, v) '.— / Vu t AVd = fv Vuef, (1.1) 

where and *3f are appropriate (infinite dimensional) function spaces. A (finite) 
Galerkin procedure consists in finding an approximation of u, U G Y 

A(U,V) = (f,V) VVGX, (1.2) 

where Y and X arc finite dimensional "counterparts" (usually subspaces, but may be 
not) of 'W and X and the bilinear form A an approximation of a. For example, when 
a = A (modulo quadrature) 3£ = W — Hj(J?) and X = Y are a space of continuous 
piecewise p-degree polynomial functions on a partition of fi, we obtain the standard 
conforming mesh-refinement (h-version) finite element method of degree p. 

The reason behind the FEM's success in such a framework is twofold: (1) the 
weak form is suitable to apply functional analytic frameworks (Lax-Milgram The- 
orem or Babuska-Brezzi-Ladyzenskaya condition, e.g.), and (2) the discrete func- 
tions need to be differentiated at most once, whence weak smoothness requirements 
on the "elements" . 

In this article, we depart from this basis by considering second order elliptic 
boundary value problems (BVP's) in nonvariational form 

find u such that A:D 2 u = f in Q and u\ dn — g, (1-3) 
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for which one may not always be successful in applying the standard FEM (with 
reference to ^2] for the notation). Indeed, the use of the standard FEM requires 
(1) the coefficient matrix A : Q — >• \R dxd to be (weakly) differentiable and (2) 
the rewriting of the second order term in divergence form, an operation which 
introduces an advection (first order) term: 

A:D 2 u = div (AVu) — (div (A)) Vm. (1.4) 

Even when coefficient matrix A is differentiable on i?, this procedure could result 
in the problem becoming advection-dominatcd and unstable for conforming FEM, 
as we demonstrate numerically using Problem (4.5 1. 

Our main motivation for studying linear elliptic BVP's in nonvariational form 
is their important role in pure and applied mathematics. An important example 
of nonvariational problems is the fully nonlinear BVP that is approximated via a 
Newton method which becomes an infinite sequence of linear nonvariational elliptic 
problems B6h08]. 

In this article, we propose and test a direct discretization of the strong form 
(1.3) that makes no special assumption on the derivative of A. The main idea, 



is an appropriate definition of a finite element Hessian given in j ]2.5| The finite 
element Hessian has been used earlier in different contexts, such as anisotropic 
mesh generation [AV021 ICSX071 |VMD+Q7j and finite element convexity |AM08j . 
The finite element Hessian is related also to the finite element (discrete) elliptic 
operator appearing in the analysis of evolution problems |Tho06j . 

The method we propose is quite straightforward, and we are surprised that it 
is not easily available in the literature. It consists in discretizing, via a Galerkin 
procedure, the BVP (1.3) directly without writing it in divergence form. 

The main difficulty of our approach is having to deal with a somewhat involved 
linear algebra problem that needs to be solved as efficiently as possible (this is 
especially important when we apply this method in the linearization of nonlinear 
elliptic BVP's). We overcame this difficulty in ^3j by combining the definition of 
it's distributional Hessian, 

(D 2 u\cl)) = -(Vu(g)V(t>) + (Vu(g)n(t)) di2 V0GC°°(I2), (1.5) 



with equation (1.3) into a system of equations that are larger, but easier to handle 
numerically, once discretized. 

It is worth noting that there are alternatives to our approach, most notably the 
standard finite difference method and its variants. The reason we are interested 
in a Galerkin procedure is the ability to use an unstructured mesh, essential for 
complicated geometries where the finite difference method leads to complicated, 
and sometimes prohibitive, modifications (especially in dimension 3 and higher), 
and the potential of dealing with adaptive methods, using available finite element 
code. Furthermore, our method has the potential to approach the iterative solution 
fully nonlinear problems where finite difference methods can become clumsy and 
demanding |KT921 ILR051 IQbeOSI ICSHH] . 

This paper focuses mainly on the algorithmic and linear algebraic aspects of the 
method and is set out as follows. In Sj2]we introduce some notation and set out the 
model problem. We then present a discretization scheme for the model problem 
using standard conforming finite elements in C°(42). In fj3] we present a linear 
algebra technique, inspired by the standard Schur complement idea, for solving 
the linear system arising from the discretization. Finally, in Sj4] we summarize 
extensive numerical experiments on model linear boundary value problems (BVPs) 
in nonvariational form and an application to quasilinear BVP in nonvariational 
form. 
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v\\l := \\v\\l hfrn = ll D " v l| 2 ( 2 - 2 ) 



2. Set up 

2.1. Notation. Let Q C R d be an open and bounded Lipschitz domain. We denote 
L 2 (i7) to be the space of square (Lebesgue) integrable functions on Q together with 
it's inner product (v,w) :— J^vw and norm ||u|| := ||v||l 2 (,q) = (^j 11 ) 1 ^ 2 - We also 
denote by (f) u the integral of a function / over the domain uj and drop the subscript 
for u = Q. 

We use the convention that the derivative Du of a function u : fl — > R is a row 
vector, while the gradient of u, Vu is the derivative's transpose, i.e., Vw = (Dw) T . 
We will make use of the slight abuse of notation, following a common practice, 
whereby the Hessian of u is denoted as D 2 u (instead of the correct VDu) and is 
represented by a d x d matrix. 

The Sobolev spaces |Cia781 IEva98j 

u k {n) : = w k 2 (n) = J <t>& L 2 (n) ■. ^ D a </» e l 2 (x?) I , (2.1) 

[ \a\<k J 

are equipped with norms and semi-norms 

2 

lH fc (r2) 

|ct|<A; 

and \v\ 2 k := = ^ ||D««|| 2 (2.3) 

\a\=k 

respectively, where a. = {ai, ad} is a multi-index, \a\ = Y2i=i a i ano - derivatives 
D a are understood in a weak sense. We pay particular attention to the cases 
fc = l,2, 

Hj(i?) := closure of C^(f2) in H 1 ^) (2.4) 
and H _1 (/2) := dual (Hj(17)) . (2.5) 

We denote by (v \ w) the action of a distribution v on the function w. If both 
v, w £ L 2 (i7) then (v\w) = (v, w). 

We consider the following problem: Find u € Hq(X?) such that 

&u = f in Q, 

u = on dQ, ^'^ 

where the data / : Q — > R is prescribed and «5f is a general linear, second order, 
uniformly elliptic partial differential operator. Let A £ L OCl ((l) dxd n Sym(IR dx ' J ), 
the space of bounded, symmetric, positive definite, d x d matrixes. 



u h-^ jSfu := A:D 2 u, 



(2-7) 



we use X.Y := trace (X T ^) to denote the Frobenius inner product between two 
matrixes. 



2.2. Discretization. Let & be a conforming triangulation of S7, namely, 3? is a 
finite family of sets such that 

(1) K € & implies K is an open simplex (segment for d = 1, triangle for d = 2, 
tetrahedron for rf = 3), 

(2) for any K, J £ £? we have that K n J is a full subsimplex (i.e., it is either 
0, a vertex, an edge, a face, or the whole of K and J) of both K and J and 

(3) \] Kesr K = Ji. 
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The shape regularity of S? is defined as 

M (^):= inf \ (2.8) 

where px is the radius of the largest ball contained inside K and hx is the diameter 
of K. We use the convention where h : S7 — > R denotes the meshsize function of 
i.e., 

h(x) := m&xtiK- (2-9) 

We introduce the finite element spaces 

V:= {$ GH 1 ^) : e P p Vife (2.10) 
V:=VnHj(i7), (2.11) 

where P fc denotes the linear space of polynomials in d variables of degree no 
higher than a positive integer k. We consider p > 1 to be fixed and denote 
by N := dimV and N = N + N 9 := dimV. Let <t> = &^) T and * = 

($!,..., 4^, $ 1 ,...,$ A r a ) T where and $ 1( $atJ forma 

basis of V, V respectively. 

Testing the model problem (2.6) with cf> € H (l7) gives 



(j^,0) = (A:D 2 U ,0) = (/»- (2.12) 



In order to discretize (2.121 with V we use an appropriate definition of a Hessian of 
a finite element function. Such a function may not admit a Hessian in the classical 
sense, so we consider it as a distribution (or generalized function) which we recall 
the definition. 

2.3. Definition (generalized Hessian). Let n : dfl — > R d be the outward pointing 
normal of fl. Given v G H (.!?) its generalized Hessian defined in the standard 
distributional sense is given by 

<D 2 w|0) = - (Vv® V0) + (Vu®n <j)) dn V</> € C°°(J?), (2.13) 

where we are using x ® y := xy J to denote the tensor product between two geo- 
metric vectors x and y. 

2.4. Theorem (finite element Hessian). For each V G V there exists a unique 
H[V] G V dxd such that 

(H[V],<f>) = (D 2 F|$) V$eV. (2.14) 



Proof . Given a finite element function FgV, Definition 2.3 implies 

(D 2 V {(/)) = - (W (8) V0) + (W ® n G C°°(ft). (2.15) 

We fix V and let 

G: C°°(I2) R d * d 

i ^ -(VV®V0) + (W®n <j)) dn . K ' 

Notice that G is a bounded linear functional on C°°(J2) in the H 1 (i7)-norm as, 

\G{<j>)\ = |(VV® V0)| + |<W ® n 0) sn | < C(d, fl) 11^ . (2.17) 

Thus, due to the density of C°°(Q) in H 1 ({2), G admits a unique extension, G. 

Let R = G be the restriction of G to V. Since G is linear and bounded on 

v 

H 1 (i7) it follows that R is linear and bounded on V in the H 1 ^) -norm. Hence by 
Riesz's Representation Theorem there exists an H[V] G \ dxd such that for each 
$ G V 

(H[V], $) := i?($) = - (W ® V*) + (W ® n , (2.18) 

which coincides with the generalized Hessian (cf. Definition 2.3) on V. □ 
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2.5. Definition (finite element Hessian). From Theorem 2.4 we define the finite 
element Hessian as follows. Let V G V then 



{H[V], *) := - (W ® V$> + (W ® n V $ e 

It follows that if is a linear operator on V. 



(2.19) 



Taking the model problem (2.12| we substitute the finite element Hessian di- 
rectly, reducing the space of test functions to V, we wish to find U G V such that 

(A:H[U],$\ = (/,$) VleV. (2.20) 

2.6. Theorem (nonvariational finite element method (NVFEM)). The nonvari- 
ational finite element solution for the model problem's discretization (2.20) is given 

(2.21) 



as U = 4> u, where u G R is the solution to the following linear system 

d d 

Du := E B a,p M -1 C ai/J u = f. 



a=l 3=1 



The components of (2.21) are given by 

B a ^ := ^<t>,A Q '' 3 * T ^ G R^ xAr , 
M := (*, <D T ) G R NxN , 



f :=(/,*> e 



3 JV 



(2.22) 
(2.23) 
(2.24) 

(2.25) 



Proof . Since H[U] G V dxd for each = 1, . . . , d , H a ^[U] = J h a ^. Then, 
testing ( 2.20[ ) with <t>, 

/>) = EE(^ ff «^]>) 

a = l 8=1 

= x;E(*' Aa,/s * Th «.^ 

= EE(*^* T )^- 

a=l /3=1 



(2.26) 



(9 



a=l 8=1 



Utilizing Definition 2.5 for each a, P = 1 ... d we can compute h Qj fl € R , noting 



U = *u, 



= - (dp<t>, d a U) + (*n fl , d a U) dn 



(2.27) 



9,3*, <9 a 4> 



Using the definition of C Q p (2.241 and M (2.23) we see for each a,/3 = 1 . . . d 

Mh aJ 3 = C a> 8U 



'a,B 



Substituting h a ,8 horn (2.28) into (2.26) we obtain the desired result. 



(2.28) 
□ 
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2.7. Example (for d = 2). For a general elliptic operator in 2-D, the formulation 
( |2.2ip takes the form 

(B 1,1 M _1 C 1 . 1 + B 2 ^ 2 M" 1 C 2! 2 + B 1 ' 2 M _1 C 1)2 + B 2,1 M _1 C 2 ,i) u = f (2.29) 



3. Solving the linear system 



3.1. Remark ((2.211 is difficult to solve). Looking at the full system setting 



D = B Q:P M~ C Q!( g multiplying out each of the matrixes and proceeding to 

solve Du = f the resulting system would not be sparse forcing the use of direct 
solvers. 



In this section we will present a method to solve formulation (2.21) in a gen- 
eral setting. This method makes use of the sparsity of the component matrixes 
B^.C^andM. 

3.2. Remark. An interesting point of note is that if the mass matrix M were 
diagonalized, by mass lumping, then for each a and /3 the matrix B a ^M _1 C ai/ 3 
would still be sparse (albeit less so than the individual matrixes B Q ^ and C a ,p). 
Hence the system can be easily solved using existing sparse methods. However mass 
lumping is only applicable to P 1 finite elements. For higher order finite elements it 
would be desirable to exploit the sparse structure of the component matrixes that 
make up the system. 



3.3. A ge neralized Schur complement. We observe the matrix D in the sys- 
tem (2.21) is a sum of Schur complements B a ' /3 M~ 1 C a ^. With that in mind we 
introduce the (d 2 + l) 2 block matrix 



M 





B 11 




M 




B 12 






M 


1 



o 
o 

o 

M 



-Ci j2 

~Cd,d-l 
— C-d,d 





(3.1) 



3.4. Lemma (generalized Schur complement). Given 

v = (hi,i, hi i2 , • • • , hd,d-i, h d:d , u) T , 
b = (0,0...,0,0,f) T , 



solving the system 



is equivalent to solving 



for u. 



d d 
a = l fl=l 

Ev = b. 



(3.2) 
(3.3) 

(3.4) 

(3.5) 



Proof . The proof is just block Gaussian elimination on E. Left-multiplying the 
first d 2 rows by M _1 yields 



(3.6) 



1 











-M- X C M " 




h M 




" " 





1 








— M _1 Ci2 




hi, 2 













1 





M - x Q, d r 




















1 


M 'Cad 




hd,d 







B 11 


gl,2 











u 




f 
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Multiplying the i-th row by the i-th entry of the (d 2 + l)-th row for i = 1, 



B 



B 



1.1 














B 1,2 • 














gc£,(i— 1 
















1.1 


B l,2 _ 


gc£,(i— 1 





-B 1 ' 1 M _1 Ci ) i 
B 1 2 M _1 Ci2 



_B d,<J-1 M _1 C ( j ) d_i 

-c2,d 



B d,d M 1 C 











" " 




hi, 2 



























u 




f 






(3.7) 



Subtracting each of the first d 2 rows from the (d 2 + l)-th row reduces the system 
into row echelon form. 



B 



B 



G 

G 



B 






d,d-l 








B 11 M~ 1 Ci 
B 1 2 M _1 Ci 



B 



d,d- 








B 



d.ii 



-B d ' d M- x C rf , d 



hii 




" " 


hi, 2 







d,d-l 







h d ,d 







u 




f 




(3 



□ 



3.5. Remark (structure of the block matrix). In fact this method for the solution 
of the system Du = f is not surprising given the discretization presented in the proof 
of Theorem |2.6| is equivalent to the following system: 



Find U e V such that 



(H[U],$) = - (VU(g> V*) + (VU® n $) ar2 V$e 
'A:H[U],$) = (f,§) V$eV. 



(3.9) 



3.6. Remark (enforcing non-trivial Dirichlet boundary values). Given addi- 
tional problem data g € H 1 ^ 2 (!7), to solve 



_Sflt =/ in Q, 
u =g on Oft, 



(3.10) 



it is not immediate how to enforce the boundary conditions. If we were solving the 
full system Du = f , we could directly enforce them into the system matrix. 

Since g G H 1 ^ 2 (I2) by an embedding it is continuous and can be approximated by 
the Lagrange interpolant with optimal order. To enforce the Dirichlet boundaries 
we introduce a further block representation 



1 


" 










. E 3 


E 




V 




b 



(3.11) 
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where E,v and b arc defined as before and Eg,vg and bg arc defined as follows 



M 





M 



















M 





... 


M 


B u 


B 1 ' 2 ... 


1 Qd.d 


v a = 


h l,U h l,2) ■ 


•> h d,d-l> h d,d> u 




ba = [0,0, 


...,0,0,g] T . 



-1,1 

-d 

-1,2 



-d 

~d,d-l 

-c d 



(3.12) 





(3.13) 
(3.14) 

Let <l>(3 = {$i, . . . , $AT a }, then the components of Eg and bg are defined as follows 
C 9 a,0 = - (^*,a Q *a T ) + a 4» e T } an e R NxNa , (3.15) 

gr9(#,£R" s , (3-16) 
where is the Lagrange node associated with $j. 



The block matrix (3.11 1 can then be trivially solved 

Ev = b E a b a . 



(3.17) 



3.7. Remark (storage issues). We will be using the generalized minimal residual 
method (GMRES) to solve this system. The GMRES, as with any iterative solver, 
only requires an algorithm to compute a matrix-vector multiplication. Hence we 
are only required to store the component matrixes B a,/3 , C Q .^ and M. 

3.8. Remark (condition number). The convergence rate of an iterative solver 
applied to a linear system Nv = g will depend on the condition number k(N), 
defined as the ratio of the maximum and minimum eigenvalues of N: 

K (N) = ^ (3.18) 

^min 

Numerically we observe the condition number of the block matrix k(E) < Ch~ 2 
(see Table [l}. 

4. Numerical applications 

In this section we study the numerical behavior of the scheme presented above. 
All our computations were carried out in Matlab® (code available on request). 

We present two linear benchmark problems, for which the solution is known. 
We take f2 to be the square S = (—1, 1) x ( — 1, 1) c R 2 and in the first two tests 
consider the operator 

1 b(x) 
b(x) a(x) 



(4.1) 



A(x) = 

varying the coefficients a(x) and b(x). 

4.1. Test problem with a nondifferentiable operator. For the first test prob 



lem we choose the operator in such a way that (1.4 1 does not hold, that is the 
components of A are non-differentiable on tl, in this case we take 

a(x) = {xlxl) 1 '* + 1 (4.2) 

b(x) = 0. (4.3) 

A visualization of the operator ( |4.2[ ) is given in Figure [2(a) | We choose our problem 
data / such that the exact solution to the problem is given by: 



u{x) = exp(— 10 \x\ 



(4.4) 
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We discretize the problem given by (4.2 1 under the algorithm set out in {2.2 nu- 



merical convergence results are shown in Figure [2] 

4.2. Test problem with convection dominated operator. The second test 
problem demonstrates the ability to overcome oscillations introduced into the stan- 
dard finite element when rewriting the operator in divergence form. Take 

a(x) = arctan (k{\x\ 2 - 1)) + 2 (4.5) 
b{x) = 0. (4.6) 

with K £ R + . Rewriting in divergence form gives 

A:D 2 u = div (AVu) — div (A) Vu. (4.7) 

The derivatives 

d a a(x) = ^ ^ (4.8) 



l + K[\x\-l^ 

can be made arbitrarily large on the unit circle by choosing K appropriately (see 



Figure 2(b)) 



We choose our problem data / such that the exact solution to the problem is 
given by: 

u(x) = sin (wxi) sin {11x2) ■ (4-9) 



We then construct the standard finite element method around (4.7), that is find 
U e V such that for each $eV 

(aVU, V$^> - (div {A) X7U, = (f, $\ . (4.10) 

If K is chosen small enough the standard finite element method converges optimally. 
If we increase the value of K oscillations become apparent in the finite element so- 
lution along the unit circle. Figure [^demonstrates the oscillations arising from this 
method compared to discretizing using the nonvariational finite element method. 

Figure [3] shows the numerical convergence rates of the nonvariational finite ele- 
ment method applied to this problem. 

4.3. Test problem choosing a solution with nonsymmetric Hessian. In this 
test we choose the operator such that b(x) is non-zero. To maintain ellipticity in this 
problem we must choose a{x) such that the trace of A dominates it's determinant. 
We choose 

a{x) = 2 (4.11) 
b(x) = {x\xl) l '\ (4.12) 

We choose the problem data such that the exact solution is given by 

{X 1 X 2 (X 2 1 -X 2 2 ) / Q 

x i+ x 2 r (4.13) 

x = 0. 

This function has a nonsymmetric Hessian at the point 0. The nontrivial Dirichlet 
boundary is dealt with using Remark 3J3 Figure [5] shows numerical results for this 
problem. 
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4.4. Test problem with quasilinear PDE in nondivergence form. The prob- 
lem under consideration in this test is the following quasidinear PDE arising from 
differential geometry: 



div 



'1 + lVul 



1 + lVitl 



(4.14) 



where y 1 + |Vu| 2 is the area element. Here we are using |Vu| 2 = DwVti. Applying 
a fixed point linearization given an initial guess u° for each n £ N we seek u n such 
that 



div 



Vu" 



/ 



'1 + |Vu 



n-l I 



(4.15) 



IVu 



n-l| 



Applying a standard finite element discretization of (4.15) yields: Given U a £ 
for each n £ N find U n £ V such that for each leV 



vu r - 



i + ivt/"- 1 ! 



,V$ = 



/ 



1 + IVJ7™- 1 ! 



(4.16) 



In fact we can work on this problem combining the two nonlinear terms. To do 
so we must first rewrite (4.14| into the form A(u, Vu):D 2 u = /. 




(4.17) 



=Au - 



DiiD 2 uVu 
l + |Vu| 2 
VuDu \ 
l + |Vu| 2 j 



:D 2 m. 



Applying a similar fixed point linearization given an initial guess u° for each n £ N 
we seek u n such that 



Vu^Du™" 1 



2„,n 



1 + IVu 



n-ll 



f 



(4.18) 



Discretizing the problem is then similar to that set out in Section 2.2 The 
component matrixes M and C a t p are problem independent, B are defined as 



B' 



ggg^ggg^ Q 
' l+ivt/"- 1 ! 2 



(4.19) 



for q^/?. 
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Figure 1. A visualization of the coefficient of the operators (4.2) 



(on the left) and (4.5 1 (on the right) 





(a) The function (x\x^} x l z + 1 over 
Q. Note the derivatives are singular 
at x\ = and X2 = 0. 



(b) The function 

arctan (5000(|a;| 2 - 1)) over Q. 
Note the derivatives are very large 
on the unit circle. 



Figure 2. Test |4.f| Errors and convergence rates for the NVFEM 
applied to a non-divergence form operator (4.2), choosing / appro- 
priately such that u(x) = exp ( — 10 |a;|). The convergence rates are 
optimal, that is for P 1 -elements (on the left) \\u — U\\ = 0(h 2 ) and 
\u - U\ 1 = 0(h). For P 2 -elements (on the right) \\u - U\\ = 0(h 3 ) 
and \u-U\, = 0(h 2 ). 




Table 1 . Test |4.1| On the condition number of E upon discretiz- 



ing problem (4.2) using P 1 hnite elements. As claimed in Remark 



dimV 


h 


K(E) 


h~' J K(E) 


16 


0.4714 


4.904 x 10 1 


10.898 


64 


0.202 


6.594 x 10 2 


26.952 


256 


0.0943 


3.665 x 10 3 


32.633 


1024 


0.0456 


1.722 x 10 4 


35.833 


4096 


0.0224 


6.894 x 10 4 


34.737 


16384 


0.0111 


3.383 x 10 5 


41.949 


65536 


0.0055 


1.337 x 10 6 


40.43 



12 



OMAR LAKKIS AND TRISTAN PRYER 



Figure 3. Test |4.2| Errors and convergence rates for the NVFEM 
applied to a non-divergence form operator ( |4.5[ ) with K = 5000, 
choosing / appropriately such that u(x) = sin (nxi) sin (tt^)- The 
convergence rates are optimal, that is for P ^elements (on the left) 
||u - U\\ = 0(h 2 ) and \u - U\ 1 = 0(h). For P 2 -elements (on the 
right) \\u - U\\ = 0(h 3 ) and \ U -U\ l = 0{h 2 ). 










— ft— X = H 1 


EOC - 2.36^* 




*S. EOC - 2. 1 82? — 




^ Ss S. EOC = 2.090? 




EOC = 3.631B*"--*. 


EOC = 2.Q453 


EOC = 3.2922*'*»n s ^ 




EOC = 3.140? 





(a) P 1 -elements 



(b) P 2 -elemcnts 



FIGURE 4. Test 14.21 On the left we present u — U Loc (k) plotted 
on a logarithmic scale as a function over u. This represents the 
maximum error of the standard FE-solution, U, to problem (4.5) 
with 16384 DOF's (h — 1/32). Notice the oscillations apparent on 
the unit circle. On the right we show ||u — f/|| L , K \ plotted on a 
logarithmic scale as a function over i?, the maximum error of the 
NVFE-solution, U, to problem (TO) with 16384 DOF's (h = 1/32). 




-1 -1 



-1 -1 



Table [2] compares the two linearizations (4.15) and (4.18). Figure [6] show asymp- 
totic numerical convergence results for NVFEM applied to ( 4.18[ ). 
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Figure 5 . Tes t|4.3[ Errors and convergence rates for the NVFEM 

on an operator (4.111, choosing / appropriately such that u(x) = 

< 2 — 2 ^ 

XlX ^2^ x 2 X2 if x 7^ 0, or u(x) = otherwise. The convergence rates 

are optimal, that is for P 1 -elements (on the left) \\u — U\\ — 0(h 2 ) 
and \u — U\ 1 — 0(h). For P 2 -elements (on the right) \\u — U\\ = 
0(h 3 ) and \u- U\ x = 0(h 2 ). 




Table 2. Test |4.4| Comparison of the fixed point linearization 
in variational form (4.151 and in nonvariational form (4.18). We 
fix / appropriately such that u(x) = sin (ttxi) sin ( ttx2)- Taking 
initial guesses U° = U° = we discretize problem (4.14) using a 



standard FEM and using the NVFEM. Denoting Ui and Ui to be 
the NVFE-solution FE-solution respectively we run both lineariza- 
tions for until a tolerance ||E/ n +i — U n \\ (resp. U n +i — U n ) < h 2 
is achieved. We compute both the stagnation point — which is the 
iteration at which the prescribed tolerance is achieved — and the to- 
tal CPU time. Notice there is significant savings in the number of 
iterations required to reach the stagnation point using the NVFEM 
over the standard FEM, however each iteration is computationally 
more costly using the NVFEM since the system is larger and more 
complicated to solve. The CPU cost for the entire algorithm is 
comparable for each fixed h. 





h 


V2/5 


V2/10 


V2/20 


v^O 


V2/80 




FEM 


Stag. Point 
CPU Time 


5 
0.50 


13 
4.02 


16 
17.51 


26 
117.58 


32 
796.58 


36 
5308.81 


NDFEM 


Stag. Point 
CPU Time 


4 
0.72 


6 
3.40 


7 

16.49 


8 

97.93 


10 
838.8 


12 
5256.84 
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Figure 6. Test |4.4| Errors and convergence rates for NVFEM 
applied to (4.141, a quasi- linear PDE under a fixed point lineariza- 
tion. We fix / appropriately such that u(x) = sin (ttxi) sin (17x2), 
taking an initial guess u° = 0. The convergence rates are opti- 
mal, that is for P 1 -elements (on the left) \\u — U\\ = 0(h 2 ) and 
|u - U\ 1 = 0(h). For P 2 -elements (on the right) ||u - U\\ = 0(h 3 ) 
and \u- U\ x = 0(h 2 ). 
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